MPSS Summer 2019
Metrum Research Group
. # A tibble: 1 x 1
. Median_CFB
. <dbl>
. 1 -14.4
mrgsolveR package for simulation from ODE-based models
C++ inside model specification formatODEPACK / DLSODA (FORTRAN)
RID, \(\eta\), \(\varepsilon\)F, ALAG, SS etc, handled under the hoodR, C++, Rcpp, boost, RcppArmadilloR is it’s natural habitatmrgsolve started as QSP modeling toolR front end to deSolveC++ interface to DLSODARcpp - vectors, matrices, functions, environments, random numbersboost - numerical tools in C++C++ code (functions, data structures, classes)SBML to mrgsolve using R bindings to libSBMLGitHub site: https://github.com/metrumresearchgroup/mrgsolve
Issues and questions: https://github.com/metrumresearchgroup/mrgsolve/issues
mrgsolve website: https://mrgsolve.github.io
User Guide: https://mrgsolve.github.io/user_guide
Vignettes: https://mrgsolve.github.io/vignettes
The first / simplest workflow.
that come from?event in this exampleplot, mutate, as_data_frame etc …%>% %>% %>% Better.
%>% ...
---------------- source: pk1.cpp ----------------
project: /Users/kyleb/Rli...gsolve/models
shared object: pk1-so-d2b54722a2a3
time: start: 0 end: 192 delta: 0.2
add: <none>
compartments: EV CENT [2]
parameters: CL V KA [3]
captures: CP [1]
omega: 0x0
sigma: 0x0
solver: atol: 1e-08 rtol: 1e-08 maxsteps: 5000
Model parameters (N=3):
name value . name value
CL 1 | V 20
KA 1 | . .
CL, V2, Q etcTHETAs)WT, EGFR etc)
Model initial conditions (N=2):
name value . name value
CENT (2) 0 | EV (1) 0
Quiz:
. [1] "/Users/kyleb/Rlibs/mrgsolve/models"
.
.
. --------------- source: effect.cpp ---------------
.
. project: /Users/kyleb/Rli...gsolve/models
. shared object: effect-so-edd44d648855
.
. time: start: 0 end: 36 delta: 0.1
. add: <none>
.
. compartments: GUT CENT PERIPH CE [4]
. parameters: VC KA K10 K12 K21 E0 EMAX EC50 KEO [9]
. captures: EFFECT CP [2]
. omega: 0x0
. sigma: 0x0
.
. solver: atol: 1e-08 rtol: 1e-08 maxsteps: 5000
. Error : the model file simple.cpp does not exist.
. [1] "Error : the model file simple.cpp does not exist.\n"
. attr(,"class")
. [1] "try-error"
. attr(,"condition")
. <simpleError: the model file simple.cpp does not exist.>
Summary:
mread to read, compile and load a modelWe use mread() to get a model loaded into R
Next, we’ll need
model %>% %>% simulate %>% take-a-look
Event object = quick / easy way to implement dose or other intervention
. Events:
. time cmt amt evid
. 1 0 1 100 1
time, evid, cmtev(...)evidevid 1, with rate==0)evid 1, with rate > 0)evid 2)
evid 3)evid 4)evid 8). time cmt amt evid
. 1 0 1 100 1
We will use a
model %>% intervention %>%
. Model: effect
. Dim: 362 x 8
. Time: 0 to 36
. ID: 1
. ID time GUT CENT PERIPH CE EFFECT CP
. 1: 1 0.0 0.00 0.000 0.0000 0.0000 157.0 0.000
. 2: 1 0.0 100.00 0.000 0.0000 0.0000 157.0 0.000
. 3: 1 0.1 91.21 8.443 0.1552 0.2225 155.7 3.460
. 4: 1 0.2 83.19 15.502 0.5817 0.8048 152.8 6.353
. 5: 1 0.3 75.88 21.354 1.2272 1.6382 149.6 8.752
. 6: 1 0.4 69.21 26.156 2.0463 2.6356 146.6 10.719
. 7: 1 0.5 63.13 30.046 3.0001 3.7278 144.1 12.314
. 8: 1 0.6 57.58 33.146 4.0553 4.8607 142.2 13.584
We need to know this to work the example; more details in a bit
Background
File name:
What would you like to “fix” in this plot?
0. [1] 0 6 12
model %>% intervention %>% simulate %>%
. # A tibble: 362 x 9
. ID time GUT CENT PERIPH CE EFFECT CP arm
. <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
. 1 1 0 0 0 0 0 157 0 1
. 2 1 0 100 0 0 0 157 0 1
. 3 1 0.1 91.2 8.44 0.155 0.222 156. 3.46 1
. 4 1 0.2 83.2 15.5 0.582 0.805 153. 6.35 1
. 5 1 0.3 75.9 21.4 1.23 1.64 150. 8.75 1
. 6 1 0.4 69.2 26.2 2.05 2.64 147. 10.7 1
. 7 1 0.5 63.1 30.0 3.00 3.73 144. 12.3 1
. 8 1 0.6 57.6 33.1 4.06 4.86 142. 13.6 1
. 9 1 0.7 52.5 35.6 5.18 5.99 141. 14.6 1
. 10 1 0.8 47.9 37.4 6.36 7.09 139. 15.3 1
. # … with 352 more rows
If we didn’t modify the output
. Model: effect
. Dim: 361 x 8
. Time: 0 to 36
. ID: 1
. ID time GUT CENT PERIPH CE EFFECT CP
. 1: 1 0.0 0 0 0 0 157 0
. 2: 1 0.1 0 0 0 0 157 0
. 3: 1 0.2 0 0 0 0 157 0
. 4: 1 0.3 0 0 0 0 157 0
. 5: 1 0.4 0 0 0 0 157 0
. 6: 1 0.5 0 0 0 0 157 0
. 7: 1 0.6 0 0 0 0 157 0
. 8: 1 0.7 0 0 0 0 157 0
mod %>% mrgsim() %>% mutate()mod %>% mrgsim() %>% filter()mod %>% mrgsim() %>% group_by()mod %>% mrgsim() %>% as_tibble()Demo: try 0.5 g Q 12 h over 0.5 hour for 1 day in normal and impaired renal function
File name:
On the fly
Persistent
mrgsim will call update for you (on the fly)
start, end, delta, addatol, rtolhmax, maxsteps, mxhnil, ixpr$OMEGA, $SIGMAtscale (rescale the output time)digits%>% %>% %>% model %>% intervention %>% %>% simulate %>% take-a-look
This is our second simulation workflow.
mod %>% ev(amt = 100, ii = 24, addl = 3) %>%idata_set( pop ) %>% mrgsim() %>% plot()
idata_set takes in individual-level data. ID CL VC KA KOUT IC50 FOO
. 1 1 1.05 47.8 0.8390 2.45 1.28 4
. 2 2 0.73 30.1 0.0684 2.51 1.84 6
. 3 3 2.82 23.8 0.1180 3.88 2.48 5
. [1] 10
. ID CL VC KA KOUT IC50 FOO
. 1 1 1.05 47.8 0.8390 2.45 1.28 4
. 2 2 0.73 30.1 0.0684 2.51 1.84 6
. 3 3 2.82 23.8 0.1180 3.88 2.48 5
.
. Model parameters (N=3):
. name value . name value
. CL 1 | V 20
. KA 1 | . .
idata?Batches of simulations or sensitivity analyses
. ID CL
. 1 1 0.50
. 2 2 0.75
. 3 3 1.00
. 4 4 1.25
. 5 5 1.50
+ mod %>%
idata_set(idata) %>%
ev(amt = 100, ii = 24, addl = 2) %>%
mrgsim(end = 120) %>%
plot(CP ~ ., logy=TRUE)Description: sensitivity analysis in a PBPK model
File name:
data_set is the dosing equivalent to idata_set. ID amt ii addl evid cmt time
. 1 1 100 24 3 1 1 0
. 2 2 300 24 3 1 1 0
. 3 3 1000 24 3 1 1 0
data_set configuration. ID amt ii addl CL evid cmt time
. 1 1 100 24 3 0.5 1 1 0
. 2 2 300 24 3 0.5 1 1 0
. 3 3 1000 24 3 0.5 1 1 0
. 4 4 100 24 3 1.0 1 1 0
. 5 5 300 24 3 1.0 1 1 0
. 6 6 1000 24 3 1.0 1 1 0
data <- mutate(data, dose = amt)
mod %>%
data_set(data) %>%
mrgsim(carry_out = "dose", end = 120) %>%
plot(log(CP)~time|factor(dose), group = ID, scales = "same")Copy from the input data to the output data
. Model: pk1
. Dim: 2886 x 4
. Time: 0 to 192
. ID: 3
. ID time dose CP
. 1: 2 0.0 300 0.000
. 2: 2 0.0 300 0.000
. 3: 2 0.2 300 2.712
. 4: 2 0.4 300 4.919
. 5: 2 0.6 300 6.712
. 6: 2 0.8 300 8.167
. 7: 2 1.0 300 9.345
. 8: 2 1.2 300 10.296
. ID amt ii addl CL evid cmt time dose
. 1 1 100 24 3 0.5 1 1 0 100
. 2 2 300 24 3 0.5 1 1 0 300
. 3 3 1000 24 3 0.5 1 1 0 1000
. Events:
. time cmt amt evid ii addl
. 1 0 1 100 1 24 3
mod + ev = ?mod + idata_set + ev = ?mod + data_set = ?. ID amt ii addl CL evid cmt time dose
. 1 1 100 24 3 0.5 1 1 0 100
. 2 2 300 24 3 0.5 1 1 0 300
. 3 3 1000 24 3 0.5 1 1 0 1000
. 4 4 100 24 3 1.0 1 1 0 100
. 5 5 300 24 3 1.0 1 1 0 300
. 6 6 1000 24 3 1.0 1 1 0 1000
Recall that data sets are just plain old data.frames in R. Feel free to code these up in any way that is convenient for you.
In our experience, the following helper functions cover many (not every) common needs for building the data sets.
expand.evas_data_setev_repev_assignexpand.evexpand.grid … gives “all combinations”. ID amt evid cmt time
. 1 1 100 1 1 0
. 2 2 100 1 1 0
. 3 3 200 1 1 0
. 4 4 200 1 1 0
as_data_setdata <- as_data_set(
ev(amt = 100, ii = 12, addl = 19, ID = 1:2),
ev(amt = 200, ii = 24, addl = 9, ID = 1:3),
ev(amt = 150, ii = 24, addl = 9, ID = 1:4)
). ID time cmt evid amt ii addl
. 1 1 0 1 1 100 12 19
. 2 2 0 1 1 100 12 19
. 3 3 0 1 1 200 24 9
. 4 4 0 1 1 200 24 9
. 5 5 0 1 1 200 24 9
. 6 6 0 1 1 150 24 9
. 7 7 0 1 1 150 24 9
. 8 8 0 1 1 150 24 9
. 9 9 0 1 1 150 24 9
ev_rep. time cmt amt evid ii addl ID
. 1 0 1 100 1 12 14 1
. 2 0 1 100 1 12 14 2
. 3 0 1 100 1 12 14 3
. 4 0 1 100 1 12 14 4
. 5 0 1 100 1 12 14 5
ev_rep. Events:
. time cmt amt evid ii addl
. 1 0 1 100 1 0 0
. 2 36 1 50 1 24 4
. time cmt amt evid ii addl ID
. 1 0 1 100 1 0 0 1
. 2 36 1 50 1 24 4 1
. 3 0 1 100 1 0 0 2
. 4 36 1 50 1 24 4 2
ev_assign. ID GROUP
. 1 1 1
. 2 2 2
. 3 3 1
. 4 4 2
. 5 5 1
. 6 6 2
e1 <- ev(amt = 100, ii = 24, addl = 9)
e2 <- mutate(e1, amt = 200)
data <- ev_assign(
l = list(e1,e2),
idata = pop,
evgroup = "GROUP",
join = TRUE
)
head(data). time cmt amt evid ii addl ID GROUP
. 1 0 1 100 1 24 9 1 1
. 2 0 1 200 1 24 9 2 2
. 3 0 1 100 1 24 9 3 1
. 4 0 1 200 1 24 9 4 2
. 5 0 1 100 1 24 9 5 1
. 6 0 1 200 1 24 9 6 2
Description: simulate GCSF PK/PD profiles for weight-based dosing
File name:
What’s going to happen?
What’s going to happen?
Combine two events
. Events:
. time cmt amt evid ii addl
. 1 0 1 200 1 0 0
. 2 24 1 100 1 24 2
What’s going to happen?
What’s going to happen?
Put two events in a sequence
. Events:
. time cmt amt evid ii addl
. 1 0 1 200 1 12 2
. 2 36 1 100 1 24 4
What is going to happen?
What is going to happen?
Wait before starting the next part of the regimen
. Events:
. time cmt amt evid ii addl
. 1 0 1 200 1 0 0
. 2 36 1 100 1 24 2
File name:
Section name:
%>% %>% %>%
mread or mread_cachemodlib("<model-name>")param(mod)init(mod)see(mod)ev(...): amt, cmt, time, ii, addl, rateDescription: Population PK of azithromycin
File name:
Description: simulate PK profiles using empirical Bayes estimates from a meropenem model fit in NONMEM
File name:
$PARAM [R] declare and initialize model parametersname = value format, or <newline>name except for words in mrgsolve:::reserved()R parsername anywhere in the modelR but read-only in the model specification file
$FIXED if appropriate$PARAM with the names in your input data sets$CMT [txt] declare model compartmentsname and number of compartments in the modelmrgsolve:::reserved()$INITExample $CMT
Example $INIT
$ODE [C++] write differential equationsdxdt_CMT = ratein - rateout;CMT for compartment amounts, parameters, or other variablesnew variables in your modelC++, which is a language$MAIN, $ODE, $TABLE, $GLOBALTo get a numeric variable, use double
Other types you might use
If in doubt, use double; it’s what you want most of the time.
$MAIN [C++] covariate model & initialsFor example
In this example
TVCL, KIN, KOUT, WT, SEX were declared in $PARAMRESPC++ expressions and functionsif(a == 2) b = 2;
if(b <= 2) {
c=3;
} else {
c=4;
}
double d = pow(base,exponent);
double e = exp(3);
double f = fabs(-4);
double g = sqrt(5);
double h = log(6);
double i = log10(7);
double j = floor(4.2);
double k = ceil(4.2);Lots of help on the web http://en.cppreference.com/w/cpp/numeric/math/tgamma
#define)POPULATION simulation$OMEGA and $SIGMA to define covariance matrices for IIV and RUV, respectively
ETA(n) and EPS(m) for realized random effects drawn from $OMEGA and $SIGMA respectively
ETAs insteadmrgsolve recognizes ID as a subject indicator for simulating a new ETA
mrgsolve allows you to write an error model as a function of EPS(m) and any other calculated value in your model
$OMEGA and $SIGMA [txt]$OMEGA and $SIGMA are allowed and each may be namedmrgsolve will enforce later on$OMEGA and $SIGMA@ delimited & must be on different line from the matrix data@block)$OMEGA and $SIGMA$OMEGA and $SIGMA are allowed; each may be namedUsers are encouraged to add labels
@ macros$OMEGA and $SIGMA@ should appear at the start of the line and the entire line is reserved for options@cor @name PKTwo forms:
@block means block=TRUE@labels ECL EVC means labels=c("ECL", "EVC")$PKMODELncmt to 1 or 2depot to TRUE for extravascular dosing compartmenttrans to change the names of required parameters$CMT to declare 1 to 3 compartments as appropriate$PKMODELTo change the bioavability of doses administered to a compartment, set F_CMT in $MAIN
To add a lagtime for doses administered to a compartment, set ALAG_CMT in $MAIN
D_CMT in $MAIN
rate = -2 in your data set or event objectR_CMT in $MAIN
rate = -1 in your data set or event objectFirst time to read
Next time to read
soloc)It can be helpful to set
and
We
. ID time GUT CENT CP
. 1 1 0 0 0 0
. 2 1 1 0 0 0
. ID time CP
. 1 1 0 0
. 2 1 1 0
. ID time GUT CENT CP
. 1 1 0 0 0 0
. 2 1 0 1 0 0
. ID time GUT CENT CP
. 1 1 0 0.0000000 0.0000000 0.00000000
. 2 1 1 0.3011942 0.6782976 0.03391488
sensFun()We use the function sensFun from the package FME (FME::sensFun)
sensFun will take the parameter values we passed and manipulate them by small amounts to determine the local sensitivitymod <- modlib("irm1", end = 72, delta = 0.1)
mod %>%
ev(amt = 100) %>%
mrgsim(end = 72, delta = 0.1) %>%
plotmyfun to capture the scenariomyfun <- function(pars, mod) {
mod %>%
param(pars) %>%
ev(amt = 100) %>%
mrgsim(end = 72) %>%
as.data.frame %>%
select(-ID)
}myfun. time EV1 CENT PERIPH RESP EV2 CP
. 1 0.0 0.00000 0.000000 0.00000000 5.000000 0 0.0000000
. 2 0.0 100.00000 0.000000 0.00000000 5.000000 0 0.0000000
. 3 0.1 90.48374 9.444382 0.04780971 4.902817 0 0.4722191
. 4 0.2 81.87308 17.851273 0.18293635 4.688450 0 0.8925637
. 5 0.3 74.08182 25.323260 0.39389717 4.426064 0 1.2661630
. 6 0.4 67.03200 31.953048 0.67040259 4.150992 0 1.5976524
. $CL
. [1] 1
.
. $V2
. [1] 20
.
. $V3
. [1] 10
.
. $KA
. [1] 1
. x var CL V2 V3 KA
. 1 0.0 CP 0.000000000 0.0000000 0.000000e+00 0.0000000
. 2 0.0 CP 0.000000000 0.0000000 0.000000e+00 0.0000000
. 3 0.1 CP -0.002544135 -0.9924278 -3.331475e-05 0.9506995
sens %>%
gather(variable,value,CL:KA) %>%
ggplot(aes(x,value,col=variable)) + geom_line(lwd=1) + theme_bw()We’ll look at Sobol’s method
./model/cipro.cppSimulate 400 mg as IV infusion over 1 hour
Plot cipro concentration over 12 hours
Simulate 400 mg as IV infusion over 1 hour
Plot cipro concentration over 12 hours
Simulate 400 mg one day of a q12h dosing
Plot over 24 hours of dosing
Explore influence of patient weight, sex, and creatinine clearance on cipro PK
idata_set, event objectFrom the Haeseker paper:
Mean (+/- SD) age was 66 (+/-17) years, the mean clearance corrected for bodyweight was 0.24 l h-1 kg^1
and the mean AUC was 49 mg l^1 h .
Table 2 shows CrCL, CL, V, and AUC by total daily ciprofloxacin dose
Let’s confirm the AUC result for each of the covariate groups in the previous step.
Use PKPDmisc::auc_partial to calculate AUC.
Profiles of unbound ciprofloxacin concentration in different tissues were predicted from the developed WB-PBPK model (Fig. 3).
Kidney and lung were predicted to have higher exposures as compared to the other organs while muscle, brain and adipose were predicted to have relatively low exposures to ciprofloxacin.
Let’s look at exposure in different organs
Use ./model/cipro_conc.cpp model
ciprofloxacin 400 mg IV over 2 hours
What is the distribution of AUC after the last dose on day 3 of q12h dosing?
Use ./model/cipro.cpp
Use ./data/cipro_post.RDS
Simulate the distribution of 24-hour AUC:MIC for
Check progress on the second day of dosing